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ABSTRACT 


Recant  Interest  in  military  VTOL/STOL  aircraft  employing  unprepared 
landing  sites  has  led  to  interest  in  the  problem  of  landing  surface  ero¬ 
sion.  Surface  erosion  Is  caused  by  the  aerodynamic  forces  on  ground 
particles  existing  within  the  flow  field  of  an  Impinging  jet.  The 
inviscid  flow  field  is  discussed  and  the  viscous  ground  boundery  layer 
is  analysed  utilising  both  theory  and  available  experimental  data.  A 
mathematical  model  of  the  process  of  entrainment  of  ground  particles  is 
constructed.  Erosion  rates  in  the  form  of  erosion  profiles  are  pre¬ 
dicted  for  selected  Jet  configuration  and  types  of  terrain.  A  criterion 
for  entrainment,  due  to  both  lift  and  drag,  was  found  end  presented  for 
selected  dlstences  from  the  jet  centerline. 
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NOTATION 


S  erosion  rat*  -  ft/aoo2 

Cjj  coefficient  or  drag 

Cf  coefficient  of  akin  friction 

Dj|  .  noaale  diameter  -  ft. 

d  particle  diameter  -  in. 

t  force  -  lb. 

Ly  lift  force  due  to  the  effect  of  the  ground  •  lb. 

Lj  lift  force  due  to  (hear  -  lb. 

m  mace  -  lb  sec2/ft. 

N  rotational  flow  parameter  ■  ^ 

H  Reynolds  number  3 

n  number  of  atripa  used  in  atrip  anelyeia 
P0  atatic  preaaura  on  the  ground  lb/ft2. 

q  ahaar  velocity  -  ft/aec. 

R  radial  diatance  from  Jet  centerline  -  ft. 
r  radius  of  particle  -  in. 

T  kinetic  energy  •  ft  lb, 

U  velocity  parallel  to  the  ground  -  ft/aec. 

V  velocity  perpendicular  to  the  ground  -  ft/aec. 
x  coordinate  parallel  to  the  ground 
y  coordinate  perpendicular  to  the  ground 

y^  vertical  diatance  at  which  ■  Uy*  in  turbulent  boundary  layer 

l  height  of  noeal*  above  the  ground  -  ft. 

^  density  -  lb  aae2/ft\ 

viacoaity  -  lb  sec/ft2. 


v 


6  angle  measured  from  center  of  particle 
OOo  vorticity 

<f  boundary  layer  thickness 

Subscripts 

crit,  value  at  boundary  layer  transition 

J  indicates  value  associated  with  jet 

L  value  at  local  point  on  the  particle 
a  property  of  air 

1  value  in  inviscid  flow 

p  value  associated  with  a  particle 

c  measurements  made  on  cylinder  in  strip  analysis 
s  static 

A  ambient 

n  measurement  made  on  a  strip 

w  value  associated  with  effect  of  wall  constraint 

S  value  associated  with  shear  effect 


1.  Introduction) 

The  increasing  In  tarts  t  in,  tad  uaa  of,  VTOL  aircraft  and  other 
vehicles  with  vertically  directed  Sans  and  jets,  par Uvular ly  in  mili¬ 
tary  oparatlona,  h«|  brought  §bput  cotjcam  for  ttaa  deterioration  of 
unprepared  or  soft  landing  aitea.  The  dovnwash  impingement  caused  by 
these  aircraft  results  in  a  radially  developed  vlseoua  flow  field  on  the 
ground,  producing  forces  adequate  to  cause  entrainment  of  the  ground 
particles.  The  reasons  for  concern  with  this  problem  are  numerous.  In 
addition  to  possible  damage  to  the  aircraft  itself  end  the  erosion  of 
the  landing  surface  there  are  the  accompanying  hacards  to  ground  per¬ 
sonnel,  damage  to  ground  support  equipment,  poor  visibility  afforded  the 
pilot,  and  the  violation  of  the  military  concept  of  concealment.  The 
present  study  was  designed  to  construct  a  mathematical  model  of  the 
ground  erosion  process  to  allow  prediction  of  the  initial  erosion  rates 
for  various  types  of  terrain  for  various  nossla  velocities  and  geome¬ 
tries.  Previous  work  [l,  14,  13,  22,  23,  25^  on  the  problem  of  erosion 
due  to  downwash  impingement  has  been  primarily  an  experimental  examina¬ 
tion  of  the  problem  with  small  scale  jets. 

The  critical  mechanism  in  the  erosion  process  is  the  interaction  of 
the  ground  boundary  layer  end  the  terrain  particles.  An  analysis  of 
this  mechanism  requires  an  understanding  of  the  viscous  mixing  problem 
of  a  finite  free  Jet,  with  the  inviscid  flow  field  extending  radially 
from  the  Jet  center  line  as  a  result  of  the  presence  of  the  ground  plane. 
An  analysis  of  the  phenomenon  which  produces  entrainment  of  the  ground 
particles  is  basically  a  study  of  the  aerodynamic  forces  which  exist  in 
the  ground  boundary  layer,  as  they  are  affected  both  by  the  decaying  in- 
viscid  flow  field  above  the  boundery  layer  and  by  the  constraint  of  the 
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ground  surface,  in  analysing  the  boundary  layer  forces  available  theory 
was  utilised  whenevar  possible,  but  sines  the  radially  developed  bound¬ 
ary  layer  runs  the  gamut  of  stagnation  point  flow  to  turbulent  flow  and 
then  complete  decay,  ttttfe  are  areas  for  which  no  theory  or  exact  solu¬ 
tion  exists.  To  fill  these  holes,  experimental  data  were  used.  Where 
exact  solutions  were  used  for  the  velocity  profiles  within  the  boundary 
layer,  verification  was  made  with  experimental  results  if  data  existed. 

The  analysis  was  restricted  to  single  jets  of  uniform  dynamic  pres¬ 
sure  loading  and  circular  cross-section  under  "no  wind"  conditions  shown 
diagrammatical ly  in  Fig.  1.  An  axi-synnetric  boundary  layer  was  assumed, 
acting  over  a  terrain  of  uniform  particle  sice  and  danslty.  Table  I  in¬ 
dicates  the  various  combinations  of  nozsle  characteristics  and  radial 
stations  that  were  analysed.  The  prediction  of  erosion  rates  is  neces¬ 
sarily  restricted  to  Incipient  motion,  since  the  secondary  effect  of 
random  collision  is  not  readily  treated  by  anelysls  and  because  the  con¬ 
stantly  changing  shape  of  the  ground  surface  la  not  accounted  for  as  a 
continuing  process.  Although  this  work  is  restricted  to  consideration 
of  uniform  Jets  only,  previous  experiments  by  Morse  [23j  Indicate  that 
the  dynamic  pressure  distribution  due  to  downwash  from  rotor  blades  and 
ducted  fans  is  similar  to  that  for  a  uniform  Jet.  The  uniformity  of 
the  Jet  or  downwash  appears  to  be  a  critical  factor  only  for  ratios  of 
jet  altitude  to  diameter,  z/Djj,  lass  than  unity. 
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2.  History  of  tha  Problem, 

The  cits sic  problem  of  in  impinging  jet  on  a  normal  surface  it  not 
new.  but  its  application  to  tha  field  of  VTOL  rircrsft,  of  either  rotor 
or  J*fc  fcyP*>  hit  batn  ganarata d  only  within  tht  past  stvral  years,  Tha 
problem  of  surface  erosion  resulting  from  impinging  downwash  first  re- 
calved  attention  from  Kuhn  Q.4}.  Kuhn's  work  was  sxpsrimental  employing 
small  seala  test  equipment  to  make  dynamic  preesura  dlatributlon  surveys. 
An  experiment  was  conductsd  by  fixing  the  nostla  dynamic  pressure  end 
relslng  or  lowering  it  over  terrain  of  known  condition.  Tasting  wes 
done  with  verious  types  of  sod,  dry  end  wet  dirt,  duet,  and  sand.  Simi¬ 
lar  experimental  tests  have  been  conducted  more  recently  by  MorBe  [23}. 

The  moat  recent  work  done  on  this  problem  is  that  of  Brady  end  Lud¬ 
wig  {26}.  This  work  consists  of  extensive  experiments  to  investigate 
the  character is tics  of  tha  impinging  uniform  jet,  including  both  the 
invlecid  flow  and  the  ground  boundary  layer.  Measurements  of  boundary 
layer  velocity  profiles  ware  compared  with  theory  and  ehowed  agreement 
to  a  limited  degree.  Some  of  theea  experimental  reeults  have  been  used 
in  support  of  the  present  work.  Vidal  [28}  has  given  an  insight  into 
tha  aerodynamic  forcea  existing  within  tha  boundary  layer. 
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3.  Th«  Flow  Field. 

In  order  to  consider  the  mechenlam  of  particle  entrainment  and  ground 
erosion,  a  thorough  knowledge  of  the  flow  field  of  the  impinging  jet  la 
re^rea.-TO 

giona.  The  first  of  these  regions  exists  regardless  of  the  presence  of 
the  ground  surface,  and  is  the  region  of  viscous  decay  immediately  exter¬ 
nal  to  the  jet  nossle.  Viscous  decay  of  a  free  Jet  is  a  classical  prob¬ 
lem  under  uniform  conditions.  This  region  is  characterised  by  the  flar¬ 
ing  out  of  the  jet  as  the  jet  boundary  mixes  viscously  with  the  stationary 
field  external  to  it.  introducing  a  ground  surface  into  the  flow  field, 
normal  to  the  jet,  does  not  invalidate  the  quantitative  solution  of  free 
Jet  decay,  but  restricts  its  extent  of  usefulness.  Experimentation  has 
shown  that  the  theory  of  viscous  decay  of  a  free  Jet  is  adequate  for 
the  impinging  Jet  problem  for  nossle  heights  of  greeter  than  about  eight 
jet  diameters.  Since  in  this  work,  the  Interest  is  in  nossle  heights 
of  Z/Dn  <  g  ,  an  additional  region  of  vlecous  mixing  with  the  station¬ 
ary  boundary  must  be  studied  and  determined,  beginning  with  the  point  at 
which  the  prese.m.«  of  the  ground  surface  has  altered  the  free  Jet  charac¬ 
teristics.  Further  study  must  then  be  made  of  the  viscous  decay  of  the 
Jet  after  impingement  on  the  ground  plane.  To  date  no  satisfactory 
theoretical  method  exists  for  dealing  with  the  entire  flow  region  in 
three  dimensions.  This  problem  may  be  partially  overcome,  however,  by 
including  the  viscous  decay  region  in  describing  the  inviscid  flow  field. 

Knowledge  of  the  inviscid  .'low  field  is  necessitated  by  the  require¬ 
ment  that  the  inviscid  velocity  exlstlrg  at  the  upper  edge  of  the  ground 
boundary  layer  be  known.  The  difficulty  in  finding  a  theoretical  solu¬ 
tion  for  this  region  lies  in  the  fact  that  although  the  boundary  conditions 
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art  wall  defined,  tha  location  of  tha  fraa  atraaalina  la  not  known.  Thraa 
dimensional  thaoratioal  aolutiona  of  thia  problem  ara  limited  to  approxi¬ 
mate  mathoda.  One  of  Chase,  by  LaClerc  £l0] ,  uaaa  an  alaetrolytie  analog 
to  datarmina  tha  ahape  of  tha  fraa  atraaalina  boundary.  Zn  thia  method 
thd  ablution  of  tha  inviscid  flow  region  waa  accomplished  by  relaxation 
tachnlquaa,  predicting  velocity  and  pressure  distribution.  Such  a 
method  of  solution  could  be  made  to  approach  tha  exact  solution  to  any 
desired  degree  of  accuracy. 

In  tha  present  work  an  axact  solution  of  tha  Inviscid  flow  field 
was  not  required.  It  was  assumed  adequate  to  use  existing  experimental 
data  which  compared  favorably  with  tha  limited  theoretical  analyses  avail¬ 
able.  In  most  instances  it  waa  found  that  tha  theoretical  solutions  were 
not  in  good  agreement  with  the  experimental  data  for  R/c»N  <  ■£f-  ■  The  approx¬ 
imate  methods  of  predicting  tha  characteristics  of  tha  inviscid  flow  re¬ 
gion  by  use  of  an  equivalent  inviscid  jet  of  greater  diameter  and  decreas¬ 
ing  dynamic  pressure,  or  by  replacing  tha  Jet  with  a  cylindrical  vortex 
sheet  of  constant  radius  extending  to  the  ground,  have  proven  unsatis¬ 
factory  whan  compared  with  experimental  data.  In  view  of  this,  available 
experimental  date  were  used  in  this  work  as  a  solution  to  tha  Inviscid 
flow  field.  Fig.  2  shows  tha  result  of  empirical  equations  fitted  to 
these  data. 

The  third  flow  region  is  chat  of  tha  ground  boundary  layer.  A  thor¬ 
ough  knowledge  of  this  region  was  most  important  to  this  work  as  the  ground 
particles  ara  predominately  immersed  in  tha  boundary  layer,  and  the  mechan¬ 
ics  of  particle  movement  originate  with  the  aerodynamic  forces  resulting 
from  the  characteristics  of  the  boundary  layer. 
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4.  The  Ground  Boundary  Layar. 

Th*  key  to  the  entire  analysis  sf  ground  erosion  is  the 
modal  used  for  tha  boundary  layar.  lound  theoretical  analysis  of  tha 
boundary  layar  was  usad  in  ehla  study  when  possible.  Experimental 
results  wars  usad  whanevar  the  theory  was  non-axis tent  or  was  in  large 
disagreement  with  these  results.  For  affective  analysis  tha  boundary 
layer  may  be  divided  into  four  regions.  The  first  is  the  stagnation 
area  boundary  layer  existing  adjacent  to  the  jet  center  line  and  extend¬ 
ing  to  ft/oN  «  .5  or  to  a  station  directly  under  tha  original  free  Jet 
boundary.  This  area  is  laminar  for  the  conditions  investigated.  The 
solution  of  this  region  utilised  the  Himanez  solution  [l8l  for  the  uniform 
impinging  jet  on  a  perpendicular  wall.  Experimental  results  were  not 
available  for  this  region,  so  that  the  validity  of  this  solution  is  not 
verified  but  was  assumed.  The  Himenec  solution  in  three  dimensional 
flow,  developed  by  Homann  [l8]  is  an  exact  solution  utilising  the  Navier- 
Stokes  equations  for  axlsyametrlc  flows.  The  solution  was  developed 
explicltyly  for  stagnation  conditions  and  the  velocity  profiles  result¬ 
ing  from  this  have  been  employed  to  R /qn  0.5.  The  nondlmenslonal 
velocity  profile  is  shown  in  Fig.  3. 

The  second  region  it  the  laminar  boundary  layar  extending  from  the 
edge  of  the  original  jet  boundary  to  the  radial  station  at  which  transi¬ 
tion  begins.  A  complete  method  of  solution  for  this  laminar  region,  by 
Smith  [l9],  is  wall  verified  by  experiment  for  4/d^  a  .5  [26]],  but 
gives  little  insight  as  to  the  limit  of  the  laminar  region  and  the  be¬ 
ginning  of  transition.  This  point  has  been  approximated  by  Brady  and 
Ludwig  [2^]  from  considerations  of  neutral  stability  of  the  laminar  bound¬ 
ary  layer.  A  Reynolds  Number  at  which  the  boundary  layer  becomes  neutrally 
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stable  eta  ba  determined  [lg^  and  may  ba  converted  to  a  Jat  nos ale 
Raynolda  Number  aa  a  function  of  QfVe>(nS')  CRiT:  ,  on  thla  baala  tha 
boundary  layar  bacomaa  unatabla  at  1.0  far  tha  jat  yalnaitiaa 

and  notala  dlaaatari  analyaad  in  thla  work.  In  view  of  thla,  and  tha 


fact  that  tha  Invlacad  valoelty  generally  raaehaa  a  paak  In  tha  vicinity 
of  r/on»  1*0  (Fig.  2),  tha  baginning  of  transition  was  takan  to  ba 
*/dn  »  1.0. 

Tha  tranaltlon  region  la  avan  mora  difficult  to  daflna.  Since 
fully  turbulent  flow  axiata  theoretically  whan  tha  praaaura  gradient  on 
tha  ground  aurfaca  la  essentially  zero ,  It  was  assumed  that  tha  transi¬ 
tion  region  extended  to  the  station  where  tha  pressure  gradient  could  ba 
takan  to  ba  negligible.  Fig.  4  shows  the  static  pressure  distribution 
near  tha  ground  takan  from  an  electrolytic  analogy  of  tha  entire  flow 
field  LlO}.  Tha  station  at  which  tha  gradient ,  f  ^  > 

was  takan  to  ba  saro  was  ^/om*  2.0.  For  tha  laminar  and  transition 
regions  of  tha  boundary  layer  experimental  results  were  used  to  form  tha 
velocity  profiles.  Thla  data,  collected  by  Brady  and  Ludwig  [26]  for  a 
limited  range  of  Jat  velocities  indicate  that  tha  non-dimensional 
velocity  profile  was  almost  completely  Independent  of  the  Jet  velocity. 
The  measurements  were  taken  at  nosale  heights  In  the  range  Z/dn«o,5 
to  4.0  at  four  locations  and  for  the  radial  range  of  */dnb.5  Co  1*33 
at  six  different  stations,  These  twenty-four  configurations  are  assumed 
to  be  fairly  representative  of  the  laminar  and  transition  regions  and 
are  shown  In  Figs.  5a  through  5h.  For  purposes  of  computation,  curves 
were  fitted  to  each  group  of  data. 

The  fourth  region  in  the  boundary  layer  la  the  turbulent  regime. 

It  was  assumed  that  this  region  begins  at  2  as  this  is  the 
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station  at  which  the  pressure  gradient  la  aaauned  nagliglbla.  A  limited 
amount  of  work  haa  bean  dona  on  theoretical  aolutlona  of  turbulent  bound* 
ary  layers.  One  of  these  that  appears  applicable  to  the  impinging  jet 
la  that  of  dlauert '.'£?}•'  This  single  aolutlon  for  velocity  profiles  in 
the  turbulent  boundary  layer,  assumed  valid  to  ^/oN*  4,  which  is  the 
radial  limit  of  the  analysis,  has  baan  verified  by  experimental  results 
at  ■=  1.33  (.26] .  The  uncertainty  as  to  the  characteristics  of  flow 

decay  limit  further  extension  of  the  analysis. 

Glauert  set  up  the  boundary  layer  equations  and  found  a  similarity 
solution  in  which  the  form  of  the  velocity  distribution  across  the  jet 
does  not  vary  along  its  length.  The  aolutlon  1s  dependent  upon  the 
assumption  of  an  effective  eddy  viscosity  as  required  to  satisfy  the  law 
due  to  Blaslua  for  flow  in  a  pipe  as  well  as  Prandtl's  hypothesis  for 
free  turbulent  flow.  The  first  solution  is  considered  to  be  valid  near 
the  ground  surface ,  and  the  second  to  be  valid  above  some  determined  dis¬ 
tance  from  the  ground  surface.  Fig.  6  shows  these  results  in  the  form 
of  a  velocity  profile,  for  the  radial  distances  and  noasle  heights  con¬ 
sidered  in  the  present  work.  Xt  was  found  that  the  shape  of.  the  velocity 
profile  was  essentially  independent  of  both  the  radial  distance  and  the 

nossle  height,  but  is  dependent  upon  the  nossla  Reynolds  Number. 

* 

The  turbulent  region  is  characterised  by  a  velocity  profile  which 
approaches  a  maximum,  U  r*oo<  >  with  increasing  y  and  then  decreases  with 
further  Increase  in  y.  As  is  shown  in  Fig.  1  it  is  expected  that  the 

* 

turbulent  boundary  layer  will  continue  to  grow  as  long  as  there  exists 
an  invlscld  upper  boundary.  Since  very  little  attention  has  been  given 
to  this  phenomenon  of  flow  decay,  and  because  of  the  nature  of  the  tur¬ 
bulent  solution,  a  fictitious  boundary  layer  thickness  is  assumed  for 
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this  region.  This  lo  do flood  to  bo  the  thloluiooo  thot  oxlata  et 


U  ■  UfrNxx  *  lt  WAS  £ouad  ellAt  the  noturo  of  thla  aaauoptlon  dooa  not 
affect  the  overall  result:  appreciably. 
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3.  Analysis  of  Aarodynamic  Forces  Within  tha  Boundary  Ltyar. 

Tha  ground  boundary  layar  la  *  non-uniform  flow  that  la  charaeter- 
laad  by  a  larga  valoolty  grad lan t  or  ahaar  layar  extending  ovar  a  large 
portion  of  tha  thickness,  Zn  addition,  tha  ground  aurfaea  or  wail  bound¬ 
ary  at  tha  baaa  of  tha  flow  forma  a  oonatralnt  on  tha  atraaallnaa  ta 
thay  attempt  to  paaa  banaath  an  obataela  baddad  In  tha  ground  aurfaea. 

Thara  ara  thraa  dlatlnct  foreaa  that  exist  within  tha  boundary 
layar.  These  ara:  1)  Drag  dua  to  tha  horlaontal  velocity  component  of 
tha  flow;  2)  lift  dua  to  tha  proximity  of  tha  wall;  and  3)  lift  dua  to 
the  velocity  gradient  in  the  flow.  There  Is  no  theory  available  which 
deala  with  this  problem  in  general.  Therefore  it  was  neceaaary  to  piece 
together  an  approximate  theory  for  each  of  the  various  forces. 

Drag  Force, 

The  drag  force  was  tha  easiest  to  account  for.  To  simplify  tha 
analysis  in  a  non-uniform  flow,  a  atrip  analysis  was  made  in  which  the 
spherical  particle  waa  divided  into  a  number  of  layers  and  It  was  assumed 
that  a  uniform  flow  acted  upon  each  layar.  Knowing  tha  velocity  of  tha 
flow  on  aach  layer,  from  the  velocity  profile,  the  force  on  each  layar 
and  tha  summation  of  horlaontal  forces  were  developed.  The  calculation 
of  the  total  drag  force  employed  tha  drag  coafflr  nt  of  cylindrical 
bodlea  measured  by  Hoerner  £9]: 

Nr,  <  IOO 

IOO  <  N(V  <  IC> + 
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N*  >  10* 

■  i 

In  addition  to  the  drag  force  tho  akin  friction  foroo  waa  calcu¬ 
lated  for  taeh  layer  and  lusnad  for  tho  partlclo.  Tho  friction  foroo,  a 
roault  of  a  nlnuto  boundary  layor  offoot  on  tho  partlelo  ltaolf,  was 
ralatlvoly  amall  but  In  aome  oaaoa  algniflcant.  Por  aphoraa  In  low  Roy* 
nolda  Number  flowo  tha  akin  friction  coefficient  can  bo  approximated  by: 

Cp  a  •  66 

Tho  tendency  of  tho  velocity  gradient  to  cauao  a  moment  on  the  partlclo 
waa  neglected. 

Lift  Due  to  Wall  Proximity. 

Theorlee  accounting  for  the  effect  of  a  wall  conatralnt  in  non* 
uniform  flow  do  not  axlot,  but  thla  effect  waa  eatlmated  ualng  the  exlat* 
lng  theory  for  uniform  flow  [12].  The  effect  of  the  ground  aurfaeo  con* 
atralnt  la  to  dlatort  the  flow  over  the  aphare,  crowding  tho  atraemllnaa, 
ceualng  a  negative  lift  force.  Tha  theory  la  baaed  on  kinetic  energy 
conaldaratlona,  maintaining  equilibrium  on  a  aphera  In  the  preaence  of  a 
well.  For  a  aphere,  In  a  flow  parallel  to  the  wall,  an  upward  force  la 
required  to  maintain  equilibrium,  Indicating  that  the  aphere  muat  be 
attracted  to  the  wall. 

Tha  method  of  calculation  of  thla  lift  force  waa  almllar  to  that  of 
drag.  The  atrip  approximation  method  waa  uaed  to  determine  an  average 
local  velocity  over  a  particular  atrip.  It  waa  then  aaaumed  that  thla 
average  local  velocity  waa  that  axlatlng  In  the  vicinity  of  the  atagnatlon 
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point,  realising  that  tha  at* gnat Ion  point  ia  ahiftad  toward  tha  ragion 
of  highar  valocitiaa,  Tha  lift  than  becomea  a  function  of  tha  aphara 
oiaa  and  tha  avaraga  local  velocity,  for  a  particle  baddad  on  tha  ground 
plana.  Tha  kinetic  energy  of  a  aphara  in  a  saving  fluid  near  an  in¬ 
finite  fixed  wall  can  be  approximated  by  a  firat  order  aolution  £l2]j 

T*  V*.  C  v  ftuf  ) 

where t 

A  •  14  I  ■**  a/a  r*/yA  ) 

B  %  '4  ^  k  C, I  ♦  a/k  f  Vy 

rn^  *  ma»a  of  fluid  displaced  by  tha  aphara. 

From  Lagranga'a  aquation  for  equilibrium 

*T/»*  -o 

If  tha  aphara  la  to  ba  maintained  in  equilibrium,  tha  raault  la  a  force 
exerted  on  the  aphara t 

F*  -  «7</+  CS77sx)  - 

F y  * 

Making  the  aubatltutlona,  the  force  exerted  on  tha  aphara  in  tha  y  or 
vertical  direction  becomes: 

Fy  -  “  bT/dy 

*  ty-t  C  AVU)  -  ('/*.  AV,1  *  14.  B  u£) 
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F7  *  AVU  -  c-  'AV^  V/c  Vy  * 


-  'A.  C4*  Vaa. m  oj*  «  ) 

Assuming  that  V,  V  <A  |, 

Fy  »  9/<s+uk*  <^KrVy4 

9 

where t 

r*X*  *.ojVarrf* 

we  gatt 

-  a//fc  V  tt  rVy* 

Assuming  that  tha  par dels  is  vary  elosa  Co  eha  ground  such  that| 

y  -  r  «  1 
or 

y  «*■  r 

tha  life  dua  Co  eha  wall  constraint  bacomas: 

Lw  ■  -a/l*to,Uk*TTr* 

For  eha  strip  analysis,  summation  yialdsi 

><-*'•  *•."<■»£  ‘•'O/n 

/ 

Life  Dua  to  velocity  oradiant. 

Thars  hava  baan  a  number  of  ehaoratioal  studios  of  eha  effaces  of 
velocity  gradients  or  shear  on  aarodynaaio  forces,  Mose  of  bhssa,  how¬ 
ever;  deal  either  with  an  unbounded  flow  or  two  dimensional  flow.  In 

\ 

competing  two  dimensional  solutions  with  three  dimensional  solutions  it 


\ 
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vaa  noted  that  the  diffaranca  la  a  tens  tn  tha  thraa  dimanaional  eolu- 
tlon  daacrlblng  tha  lateral  component  o£  flow.  For  tha  axl-symmetrical 
caoa  dealt  with  here*  thla  component  vaa  aaaumed  negligible.  A  aolution 
for  non-Mnifora  flow  on  cylinders  by  Murray  and  Ml to he l l  [l3]  vaa  em¬ 
ployed  in  conjunction  with  a  atrip  analyaia. 

In  tha  atrip  analyaia,  the  particle  waa  divided  into  layera  of  cir¬ 
cular  cylinders  situated  such  that  their  centers  formed  a  Line  perpendicu¬ 
lar  to  the  flow.  Beginning  with  the  final  result  of  Murray  and  Mitchell, 
the  dimena Ion laae^ 'shear  velocity  on  each  cylinder  waa  found  to  bet 


-  yN  s\^\rs  Ume)] 
K'*n  CO  coi  Une) 

n-1  KinO)  n 

C-0n  +  l  feftnV)  <S> 


n«o 


where  t 


an  CO  “  L  l^ar\  CoO  ]  | 

rv0  “  l 


r*0  P  1  CP  ♦  v) 


{.  lo^'4-  &L  Vvre-i)  4  t;n4f4|)]} 
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and 


+  C-D0*' 


f 


n  t  An 


r»o  rj  ^n  +  'O)! 


f'C.O  a  Euler’a  Constant 
■  0.57721367 


Tha  dimensional  staaar  velocity  la  given  byj 

%'  U l 

which  la  of  tha  form  (whan  squared): 

<***  =  uk*  Ca*  +  c*  +■  */ag>  +  +  *fcc) 

Where : 


F\k  *  s*nle  [,  cosv-i7  C  sin©  -  Vn  co fV>  <.sm  © )• 
sinv^c^sm©)  +  */n *"•  si^Nn*  C  s ©)] 

®*  -  t-1  ^  ^]Tf ^ 

c, .  4CS  tr  a«g>KV,„«],i|..nWW).]‘ 

*AS*  COS  An©]  • 

r'H|  K^n^l)  J 

sin  ©  QcosK  c.  a  me)  -  s ir%x  i  v><©)] 


z* Q  *  4 L ^ 0Cr  * 1  ^  K g ^ t CO  Sm^A^i)©]  • 

S  in®  [cosK  C.sm©')-  Vn  S\nV\  C^sin  ©>  3 
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^QC  *  s/m[5-  CriY'iE&nJsi}  cf)  cosine"]  • 
n«i  CO  An  1 

^SS^K^<a> 4,!'^’i>.# I .  ■ 

4Ju  is  tli*  local  velocity  at  the  center  of  each  cylinder,  and  I,  %  and 
K*  are  Baeeel  functions  of  the  first  and  second  kind,  end  the  derivative 
of  the  second  kind  respectively.  N  is  a  flow  parameter  dependent  upon 
local  velocity,  radius  of  the  cylinder  and  the  vorticlty  at  the  stagna¬ 
tion  point,  end  is  defined  as:  N*  U^/u^  .  The  determination  of 

the  flow  parameter  was  approximated  since  the  stagnation  point  on  each 
cylinder  could  not  be  known  in  advance  of  the  calculations.  In  order  to 
keep  the  shear  velocity  in  the  direction  it  was  known  to  flow,  a  lower 
bound  of  N  »  I  could  be  established.  If  Ms  and  U  are  evaluated  at 
the  geometric  centers  of  the  cylinders,  experimentation  has  Indicated  an 
upper  bound  of  N  m  "3  £13].  After  analysing  the  computed  results  of 
shear  velocity  with  values  of  N  ranging  from  1.1  to  10,  a  value  of 
N  ■  1.5  was  accepted  as  being  representative  of  tha  physical  model. 

Having  found  the  shear  velocity  the  lifting  force  due  to  shear 
affect  was  calculated  for  each  cylinder  end  sumned  over  the  particle. 

The  resulting  expression  1st 

L$  "  ^  ’  br<-  -  ^n)  sin©  c/a 
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Where: 


r<_  m  radiua  of  cylindrical  atrip* 
b  ■»  width  of  atrip 
T\  ■  number  of  atrip* 

Noting  that) 


Sin  a-  c/» 


*  O 


Wa  have: 

5  R/k  UK*  ^  ♦  fV-*  0“  4  XftQ  +  2AC.  +  *C>Q). 
2mdo/« 

In  evaluating  the  integral  we  have: 


Since  the  laat  two  integral*  are  identically  aero.  Further  we  may  write: 


r*rr  rtfT 

J  s>n©f/»  »  /  s m*1  flp  cosK*’  Csmtt)  c^e 

0  / o 

~  f*  S'^a©  cosK  C.sm®V  - 

..  r9.rr 

*  Jo  sinV^s*^. , 
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rnr  mt 

JQ  Alsim*ec/©  ■  ~Z^NjL  Sh*®  toiVi 

SinV>  C.  *m*)  c/a 

Sine*  the  fleet  end  last  Integrals  ere  Identically  aero. 

Due  to  the  complexity  of  the  above  expressions  and  the  difficulty 
of  integration,  it  was  desired  to  cheek  this  method  with  another  method 
t.29]  of  a  much  mure  approximate  nature  for  calculating  the  lift  due  to 

shear. 

The  alternate  method  of  calculating  lift  assumes  that  the  previously 
discussed  lift  due  to  the  ground  surface  constraint  is  negligible  and 
that  the  lift  forces  on  the  particle  can  be  estimated  by  using  the  solu¬ 
tion  of  Hall  8  for  a  sphere  in  a  two  dimensional  uniform  shear  flow. 

In  addition,  the  approximate  method  assumes  that  the  sphere  rests  on  a 
bad  of  spheres  of  the  same  diameter  and  that  the  pressure  on  ths  bed  is 
the  ambient  stream  static  pressure.  The  lift  coefficient  then  becomes 
simply  a  function  of  the  local  slop*  of  the  velocity  profile,  , 

ths  local  velocity  on  the  sphere,  and  the  radius  of  the  sphere.  In  view 
of  ths  many  restrictions  placed  on  this  solution  it  was  used  only  as  an 
order  of  magnitude  check  on  ths  other  solutions.  Fig.  7  shows  a  compari¬ 
son  of  the  two  solutions. 
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6.  Partial*  Entrainment  and  Inclplant  Broalon  Rata*. 

Only  a  few  of  th*  poaaibl*  typca  of  tarrain  to  which  a  mathematical 
analyala  of  th*  aroalon  rates  night  b*  practically  applied,  w*r*  lnvesti- 


gat*d.  tlnlfom  partial**  ranging  in  aia*  fvon  ,003  in,  to  ,125  in.  in 
d tan* tar  war*  invaatigatad.  Th*  email** t  ait*a  r*pr*a*nt  looa*  dirt 
and  th*  larger  diameter*  rapraaant  aand  or  aondy  gravel,  with  .125  in. 
being  repreaentativ*  of  email  gravel  type  rock.  It  waa  aeaumed  that  all 
terrain  waa  without  any  moiatur*  content.  The  denaltiaa  of  th*  particle* 
ranged  from  75  lb/ft3  for  looae  dirt  to  125  lb/ft3  for  gravel.  Table  H 
liata  th*  varioua  terrain  particlea  analysed. 

There  are  three  possibilities  for  entrainment  of  the  ground  particle. 
The  firat  ia  that  it  may  have  enough  lift  force  on  it  to  be  lifted  directly 
from  the  ground  aurfac*  into  th*  free  atream.  The  aacond  la  that  after 
rolling  motion  la  atarted  due  to  th*  drag  force,  it  may  be  able  to  gain 
th*  additional  lift  required  for  entrainment  in  th*  free  atream.  Th* 
third  la  that  with  many  particlea  rolling  at  different  velocities  th* 
random  collision  of  particlea  that  follows  will  bounce  some  of  them  into 
the  invlscld  atream.  Only  incipient  motion  la  within  th*  acopa  of  this 
work,  hence  it  was  assumed  that  any  particle  which  was  moved  contributed 
to  th*  eroalon  ret*. 

The  net  drag  force  and  the  net  lift  force,  including  both  mechanisms 
of  lift,  were  computed  for  each  of  th*  particlea,  at  each  of  thirteen 
radial  stations  on  the  ground  aurfac*  for  each  of  th*  nossl*  configura¬ 
tions  studied.  To  calculate  the  Initial  net  drag  force,  the  previously 
mentioned  ground  surface  model  of  nested  particlea  was  assumed  using  a 
coefficient  of  static  friction  of 
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7.  Results  and  Conclusions. 

Tha  partial*  sis*  crltsrla  for  lifting  antrainsiant  and  drag  entraln- 
mant  ara  shewn  la  Figs.  8,  9,  and  10  in  which  tha  particle  slsa  is  plotted 
versus  a  load  parameter  at  various  stations  with  Jet  velocity  and  jet 
height  as  parameters.  The  load  parameter  is  defined  as  the  particle  den¬ 
sity  divided  by  the  dynamic  pressure  in  the  fra*  stream.  Tha  lift  en¬ 
trainment  criterion  predicts  that  an  optimum  particle  sit*,  tha  particle 
else  most  prone  to  entrainment,  exists  for  all  radial  stations  and  that 
the  critical  density  increases  with  radial  distance  to  a  certain  . 

The  optimum  particle  sle*  for  lift  entrainment  occurs  because  particles 
smaller  than  the  optimum  are  affected  by  much  smaller  velocities,  and 
particles  larger  than  optimum  are  affected  more  by  a  flow  region  where 
the  shear  gradient,  ,  approaches  saro. 

The  drag  criterion  indicates  that  extremely  small  particles  of  0.02 
in.  diameter  or  less,  up  to  very  high  densities  will  roll  on  the  ground 
(see  Fig.  8).  As  tha  particle  slsa  Is  increased  the  possibility  of 
rolling  decreases  rapidly  up  to  a  particle  diameter  of  0.0?  In.  at  which 
else  the  tendency  to  roll  increases  again,  since  the  larger  velocities 
in  the  boundary  layer  begin  to  act  on  the  particle,  until  an  optimum 
diameter  of  approximately  0.06  In.  is  reached.  For  radial  distances 
less  than  ~  *li  velocities  in  the  boundary  layer  are  not  suffi¬ 

cient  to  move  large  particles  of  high  density.  The  range  of  optimum  else 
particle  for  drag  antralnmant  of  approximately  .03  to  .09  is  consistent 
for  all  radial  stations.  It  is  of  particular  significance  that  the  rang* 
of  optimum  particle  diameter  for  both  lift  and  drag  is  very  nearly  equal 
to  the  boundary  layer  thickness  at  the  particular  radial  station  being 
considered. 
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The  incipient  eroalon  retea  for  aelected  denaltlaa ,  perticle  a lie a , 
iud  uiiali  »l»  fiuOwu  iu  rig*.  11  And  12.  IliiSi  mXm  {3UQ« 

by  combining  «a  vectors  tht  incipient  horlsontal  end  vertical  accelera¬ 
tions  of  the  particles.  This  was  aaauaad  to  be  proportional  to  the 
Incipient  erosion  rate.  The  erosion  rate,  normalised  by  the  maximum 
rate,  is  plotted  versus  the  nondlmensional  radial  distance.  For  all  con¬ 
figurations  and  particles  the  erosion  profile  is  maximum  in  the  area 
.5  <  <2.0.  The  maximum  erosion  for  larger  particles  occurs  near 

1.  and  then  drops  off  sharply  with  increased  radial  distance. 

As  the  particles  become  smaller  a  maximum  eroalon  occurs  closer  to  the 
perimeter  of  the  jet  nossle  and  a  second  maximum  begins  to  occur  at 
^  l»3»  This  occurs  after  the  transition  region  has  begun  to 
form  and  where  the  onset  of  turbulence  can  affect  the  particles.  Fig.  13 
shows  the  incipient  rate  of  lift  entrainment  by  itself  and  indicates  that, 
in  general,  the  shape  of  the  erosion  profile  can  be  attributed  to  thd 
lift  forces. 

The  results  shown  in  the  graphs  represent  only  a  few  configurations 
selected  from  the  large  amount  of  data  collected.  In  viewing  all  of  the 
results  and  comparing  the  many  combinations  of  configurations,  it  was 
found  that  the  variation  of  ground  erosion  rates  was  principally  a  func¬ 
tion  of  the  dimensionless  radial  distance  from  the  jet  centerline  and 
the  else  of  the  particle  for  a  given  nossle  height,  .  The  ero¬ 

sion  rate  was  nearly  proportional  to  the  particle  density  except  near  the 
stagnation  point,  where  the  erosion  rates  were  very  small. 

From  these  results  it  may  be  concluded  that  for  ground  particles  of 
less  than  .04  in.  diameter  the  incipient  erosion  rate  has  two  maximum 
points.  One  in  the  vicinity  of  the  perimeter  of  the  Jet  nossle,  and  the 
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second  at  a  radial  distance  approximately  two  dlaaiatara  from  tha  stagna¬ 
tion  point.  For  vary  small  par tic las  and  low  jat  velocities  tha  second 
maximum  is  tha  most  significant.  In  all  other  easeo  where  two  maximum 
points  exist,  the  f  ira  t  is  themost  significant.  For  ground  particles 
whose  diameters  are  greater  than  .04  in.  a  single  maximum  point  axiste 
near  tha  perimeter  of  the  jet  noasle.  As  the  particle  alee  increases 
tha  influence  of  the  turbulent  region  decreases. 

Xt  was  concluded  that  particles  of  approximately  .02  in.  diameter 
or  less  would  entrain  due  to  drag  even  at  vary  high  danaltias  or  vary 
low  velocities.  For  particles  larger  than  approximately  .02  in.,  an 
optimum  slae  of  approximately  .06  in.  exists  up  to  radial  distances  of 
four  noeala  diameters  for  drag  entrainment.  An  optimum  particle  slse 
moat  susceptible  to  entrainment  due  to  lift  axiste  in  tha  range  of  .07 
to  .12  in.  diameter  for  radial  distances  up  to  four  nosala  diameters. 

Also,  it  may  be  concluded  that  varying  the  nossle  diameter  had  very 
little  affect  on  tha  dimensionless  erosion  rate;  hence  in  the  final  analy¬ 
sis  the  noeala  diameter  was  not  considered  as  a  parameter.  Moreover, 
while  the  actual  erosion  rates  were  approximately  linearly  dependent  upon 
the  Jet  velocity,  tha  dimensionless  erosion  profile  appeared  to  be  inde¬ 
pendent  of  Jet  velocity,  with  tha  exception  of  tha  order  of  importance 
of  the  maximum  points. 
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CRITICAL  PARTICLE  SIZE  FOR  DRAG  ENTRAINMEHT 
Vj=  65  ft/sec,  Z/Djj  =  1,  fy'Dy  =»  .1|_ 


CRITICAL  PARTICLE  SIZE  FOR  LIFT  ENTRAINMENT 
Yj  =  65  ft/3CC,  Z/Dh—1,  iv'dm  =  .o5 


CRITICAL  PARTICLE  SIZE  FOR  LIFT  ENTRAIENENT 
Vj  =65  ft/sec,  Z/DH  =  1 ,  r/Djj=*.2 


u.  .6  .6  1.0  1.2 

6^/^Uj2  —  l/ft  xl0-3 
FIGURE  9d 

CRITICAL  PARTICLE  SIZE  FOR  LIFT  ENTRAINMENT 
Vj=65  ft/sec,  Z/Djj  —  1,  R/Dr  -  *3 


CRITICAL  PARTICLE  SIZE  FOR  LIFT  ENTRAINMENT 
Vj=65  ft/sec,  Z/DN=1,  R/Dn=.5 


CRITICAL  PARTICLE  SIZE  FOR  LIFT  EKTRAI NMEHT 


CRITICAL  PARTICLE  SIZE  FOR  LIFT  ENTRAINMENT 
Vj  =  65  ft/sec,  Z/Dtf  =  1,  R/Dn=1i-.0 


FIGURE  10a 

CRITICAL  PARTICLE  SIZE  FOR  DRAG  ENTRAISHEHT 
Yj  —65  ft/sac,  Z/Dfl-ii-.O,  R/Dk  -  .5 
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CRITICAL  PARTICLE  SIZE  FOR 
Vj-65  ft/sec,  Z/Djj  =  lj.,0 


*  «  asjpj'.'P”  ' 


INCIPIENT  EROSION  RATS 

PARTICLE  DIAH.=  .007  In,  Vj~65ft /see,  2./%  = 


INCIPIENT  EROSION  RATE 

PARTICLE  DIAK. ~  .01  in,  Vj=65  ft/sec,  Z/Dji 


0*1 


INCIPIENT  EROSION  RATE 

PARTICLE  DIAH.  =  .02  In,  Vj=  65  Tt/sec,  Z/%^ 


UiCIPIEKT  EROSION  RATE 

PARTICLE  DIAM.-  -Ol^  in,  Vj=6£  ft/sec,  Z/Ds  = 


1 


INCIPIENT  EROSION  RATE 

PARTICLE  DIAM .  —  .003  In,  .  Vj  =  £00  ffc/sec,  Z/Dh  = 


0*1 


INCIPIENT  EROSION  RATE 

PARTICLE  DIAH.= .007  in,  Vj  =  500  ft/sec,  Z/Djj  = 


1.0 


INCIPIENT  EROSION  RATE 

PARTICLE  DIAH.  —  .01  In,  Vj  =  500  ft/sec,  Z/D»|  = 


INCIPIENT  EROSION  RATE 

PARTICLE  DIAM.  »  .02  la,  Vj- 500  ft/sec,  Z/Dh  = 


1.0 


INCIPIENT  EROSION  RATE  FOR  LIFT  EHTRAIKKEST 
PARTICLE  DIAM.  =  .007  In,  Vj=6£  ft/sec,  Z/Du^.$ 


1.0 


INCIPIENT  EROSION  RATE  FOR  LIFT  FffTRAIHKSK' 
PARTICLE  DIAM.=  .0\  In,  Vj  =  6$  Pt/sec,  Z/%  = 


1HCIPXKHT  EROSIOH  RATE  FOR  LIFT  ENTRAIHMESI 
PARTICLE  DIAH.  =  .02  In,  Vj=65  ft/sec,  2/Djj  = 


INCIPIENT  EROSION  RATE  FOR  LIFT  ENTRAINMENT 
PARTICLE  DIAM.  — .Ol|.  In,  Vj=65  ft/sec,  Z/l)Ht= 


FIGURE  l!j.f 

1  INCIPIENT  EROSION  RATS  FOR  LIFT  ENTRAINMENT 
PARTICLE  DIAM.  =  .06  in,  Vj=  65  ft/sec,  Z/Djj^.5 


FIGURE  li|g 

INCIPIENT  EROSIOH  RATE  FOR  LIFT  ENTRAINMENT 
PARTICLE  DIAK.=  .125  in,  Yj=  65  ft/sec,  Z/Dg=*5 


FIGURE  15 

LAMINAR  BOUNDARY  LAYER  THICKNESS 
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APPENDIX 


FORTRAN  PROGRAMS 


MPT  FORCES: 

Input. 

Program  is  generated  from  the  polynomials  describing 

the  velooity  profiles  of  tha  boundary  layer*  for 

various  radial  stations. 

MM  Number  of  terms  ol'  tha  polynomial 

PCOZ  Coaffloionts  of  tha  polynomial 

QSH  Nondlmensional  shear  velooity 

VJ  Jet  velooity 

DN  Nozzle  diameter 

RG  Radial  dlstanoa  from  jet  centerline 

ZO  Nozzle  height  above  ground 

YHVX  Dlstanoa  to  have  maximum  velooity  in  turbulent 
boundary  layer 

VMX  Maximum  velooity  in  turbulent  boundary  layer 

Intermediate  Results. and  Output 

DP  Particle  diameter 

VI  Inviaoid  velooity 

YND  Nondlmenaional  dlstanoa  above  ground 

XND  Nondlmenaional  roots  of  velooity  profile 

polynomial 

DNS  Partlole  density 

DVX  Derivative  of  polynomial 

SDVX  Slope  of  velooity  profile  •• 

DUDY  Slope  of  velooity  profile  - 

FLSC  Shear  lift  {oyllnder  method) 

SLFT  Shear  lift  (sphere  method) 

PLPT  Lift  due  to  wall  oonetraint 

PLAC  Aooeleratlon  of  partlole  in  vertioal  direotlon 

Subroutine  Root  finds  the  roots  of  the  velooity  profile 
polynomial 
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SUBROUTINE  L I  FT  C  MM,PCOZ»  QSH, V J, DN, RG , ZG, NCC , YHVX, VMX > 

DIMENSION  DP ( 10  ) «  RP<  10 ) «  YQU  ( 40 ) ,  YGBUO  > ,  YL(  70  > .  YQ!  7Q  > , 

1AR<7o>,VL<70>,PCOZ(3o>#XL<4o>,FX<30>,NLY!5  ,DVX<30># 

2  DGB(5)iPCOX(10)<R(100),YnD(120),XND(120)#YNDD(50),XNDD(50). 

3  VLO(50) , Y0DD(20)  V 

DPU>«.Q6.  f 

DP(|>*»64 

DPl3)»,o2 
DP( 4 ) » ,  01 
DP(5)«,007 
DP! 6) » . 125 
DP( 7 ) #  .  003 
RGN«RG/DN 
ZGN»ZG/DN 
PRINT  611*  DN 

611  FORMAT!///  22HN0ZZLE  DIAMETER  (FT)  ■  F5.2  ) 

PRINT  612.  VJ 

612  FORMAT!  26HN0ZZLE  VELOCITY  (FT/SEC)  a  F7.1  ) 

PRINT  4Q0.RGN 

400  FORMAT (26HN0NDIMEN  RADIAL  DISTANCE  «  F2Q.10) 

PRINT  401,  ZGN 

401  FORMAT (2BHN0NDIMEN  HEIGHTH  OF  NOZZLE  *  F2Q.10) 

DO  115  I  1=1.7 

PRINT  403,  DP ( 1  I ) 

403  F0RMAT(//23hpaRTICLE  DI AMETER C I N .  )  *  F2o • 10  )  • 

'RP(II>«DP(lI)/2. 

GNU«3. 7456-7/2. 3769E-3 
IF(NCC-25)3o2,3o3,3q2 
303  VI«VMX 
GO  TO  700 

302  IF(ZG/DN-. 5)640, 640, 661 
661  KZQiZG/DN 

GO  T0(640,641,640,643) ,KZQ 

640  IF(RGN-. 35)60,60,61 
61  IFCRGN-. 75)18,16,19 

19  IF(RQN-1,2)20,20,22 

60  VI«VJ*SORTF  < ABSF ( 1 . B8*RGN**2- . 225*RGN )  ) 

GO  TO  700 

18  VI«VJ*S0RTF<-3,2*RGN**2+5,13*rqn-1.26) 

GO  TO  700 

20  VI»VJ*SQRTF(-2,68*RGN**2+5,48*RGN-1.82) 

GO  TO  700 

22  VIaVJ«SGRTF<2.25*RQN**2»9.45*RQN+8.96> 

GO  TO  700 

641  IF<RGN-,35)6q,6o,644 

644  IF(RGN-1. 0  )645, 6.45, 647  ' 

645  VI«VJ*SQHTF<-.0985*RGN**2*1.13*RGN-.241> 

GO  TO  700 

647  VI»VJ*SORTF< . 077#rgn**2-,64i#rqn*i,354> 

GO  TO  700 

643  IFCRGN-. 40)62,62,654 

654  IF(RGN-1. 2)655, 655, 657 

655  Vl«Vj*SQRTF!-.5*RGN**2+1.2*RGN-.3> 
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50  TO  700 

657  V  I*VJ*SORTf<  .  t)4o4*RGN**2-  ,  3l8*RGN*  .  744  ) 

GO  TO  700 

62  VI*VJ*SQHTF  (,625*RGN**2*.o25*Rqn> 

700  ANCLU01. 

016  TH«DP< 1 1 )/ANCL 
NNtANCi. 

NA* ( NN+1 ) 72 

B  (  1  )  a  B  B  (  J  T  ) 

DO  1  K*2,NA 
,  AA«K-l 

xiji-aa*th  . 

1  R(K>*SQfnF<RP<IJ>**2-XLC**2> 

100  If(NCC-26>901*900,900 
901  If <NCC-25)?04. 709*900 

7 0*5  YN0<  l)«Rtl>/YHVX _ 

GO  TO  10 ^ ‘ 

704  YNDC  1 >*  <R<  1 ) / < 12 . *DN > > *SQRTF ( V J*DN7QNU ) 

GO  TO  104 
900  AAA*  V I /R(3 

YND  (  1  )*H<1  >*SQRTF  <2  .  *A AA/GNIJ  >  712 , 

L 0 4  X ND  ( 1 )  * 0  .  I) 

IF(NCC-25>50S.504,508 
504  I F( YND  <  1  )-. 15)500  *500  *501 
500  MM*8 

DO  511  I  2*1 »  MM 
511  PC0X<U)*PC02<IZ) 

GO  TO  503 
>01  MM*  4 

DO  5o2  1*1, MM 
N*  l  *8 

502  PCUXII  )*HCOZ<  N) 

506  DO  507  KW*1 »  MM 
6X1*MM-KW 

FX(KW)«PCOX(KW)*YND(  1)**£X1 

507  XND<  1>*XND<  1  > ♦FX  ( K W > 

GO  TO  673 

508  DO  71  KZ«1.MM 
EX1*MM  -KZ*1 

fX(KZ>  «PC0Z  <KZ)*YND<l)**eXl 
71  XND<l>«XND(l)+fX<KZ) 

GO  TO  673 

105  CALL  HOQI'U  , MM, PC0Z.  XND*  YND* NCC > 

573  If (XND(l)-l.O >31*31.32 
32  XND< 1 ) * • 9999 
31  VLll>«Vl*XND(l> 

706  A R T ■ 0  « 

DO  43  JZ«1*NA 
If  (JZ-D44, 44,45 

44  ARE*2**R(JZ >*3.141 59 *TH 
GO  TO  43 

45  ARE*4. *R< J2 >*3,14l59*TH 
43  ART«ART*ARb 

SPA«RP< l l >**2*4, *3.14159 
ADS*Sr  A«AR1 
ADS*  ADS/SPA 
AC0R*1 . *  ADS 
' FL5*  0 • 

DO  39  IZ*1*Na 
lfUZ*l>  40  *40.41 

40  FL*-TH/2.*2.37696-3*RUZ>*VL(l>**2*0SH/144. 
GO  TO  39 
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P  L.* 


■2*376yE-3*fiiJZ)*VL(i)**2*uSN/l44. 


39  FLS*FLS*FL 
FLSC*FLS* ACOR 

APPROXIMATE  shear  lift  OF  particle  BY  SPHERE  METHOD 
130  YGD0( 1 ) *•  H <  1 ) 

IF<NCC-25>132,77,129 
77  YNDD(1 ) ■  YQDO<  A)  /YHV,X 

GO  TO  133  "~J . . 

129  Y.NDDIl5  =  YCDD(i)*SSRTF  *  2 . ■ A A A/GNU  5 /!«  . 

90  TO  133 

1*®  VNpOU  ) « < VlAfi <  1 ) / (  I* . *DN ) ) PlSRmV : J*ON/«NU > . - .  ..... 

-f33  SBv^rg jo .  . . 

DO  74  KD"1<MM 
MX2*MM-KD 
EX1*MM-K.D*1 
JFCNCC-25)75, 79,75 
79  PCOZ(KD)aPCOX(KO) 

75  IF<KD-MM)180, 181,180 

180  DVX(KD)«EXl*PC0Z(KD)*YN0D(i)**MX2 
GO  TO  74 

181  DVX ( KD ) »PCQZ ( KD  ) 

74  sdvx«sdvx*dvx<kd> 

IF(S0VX>182, 183,183 

182  IF(NCC-25)184,183,184 
184  SOVXsO.O 

IF(NCC“25)l9Q*183*192 

183  DUDY»S0VX*VMX/YHVX*12. 

GO  TO  193 

192  DUDY«SDVX*SQRTF<2.*AAA/GNU)*VI 
GO  TO  193 

190  DUDYbSDVX*SQRTF(VJ*DN/GNU>*VI/DN 

193  20N«  R(1)/VL<1)*DUDY/12, 

CLF-0. 345*0 .557*ZON*0,879*ZON**2 

SLFT«CLF*2, 3769E-3*VL(1 J**2/2.*3,14159*R(1>**2/144, 

PRINT  800 

800  FORMAT (14HN0NDJMEN  SHEAR, 5X , 17HNONDI MEN  VELOC I T Y , 5X. 17HUI FT  DUE  TO 
1  SHEAR, 5X.17HUFT  DUE  TO  SHEAR , $X .  19HAREA  DIFFERENCE  0/0, 3x,  13HPR0 
2FILE  SLOPE) 

PRINT  801 

801  FORMAT  ( 3X.8HVELOC  I  TY,UX,HH0N  C  YL  I NDER ,  12X,  10H<  CYL I NOBR ) ,  13X, 

1  8H(SPHERE)/) 

PRINT  802 < QSH , XNQ  ( 1 ) , FLSC , SLFT , ADS, SDVX 

802  FORMAT ( 6E20 • 1 0 ) 

GO  TO  52 

5  PRINT  200 
200  FORMAT (/^HfcHROR/) 

52  NNT «59 
UN»NNT 

603  DGB ( 1 > - 0 . 

YH*0  • 

YLM«NNT-1 
YINC*DP<  1  D/YLM 
DO  169  J»1,NNT 
|F(J-i>666, 111,112 

111  YG( J ) "0 . 

GO  TO  169 

112  YG  ( J ) * YH* Y I NC 
YH*YG( J) 

169  CONTINUE 

DO  113  K«2, NNT 
■  YND( 1 ) *0 . 

| F  <  NCG-26  )  116 , 902, 90?__ 

914. 


1 


£  i 


902  YND<K)«YG<k)*SQRTF<2,*AAA/GNU>/12. 

GO  TO  113 

116  |F(NCC-25)30C*30i*902 
301  YND<!C>9Y9(K)/YHVX 
GO  TO  113 

-  i M  YQ(K  > / 112  i«DN ) .) ,6lflHTf  ( VUADN/GNU)  . . . 

113  ffONTlNu* -  -  -  . . . - . -  - 

671  DO  678  NX»i*NNT 
XND ( NX  )  ■ 0  . 

IF<NCC-25)708*720,708 

720  IFIYNDINX)-. 15)721  *721,701 

721  MM«8 

DO  811  I 2»1*  MM 
811  PCOXUZ)  «PCOZ<  IZ) 

GO  TO  708 
701.  MM«4 

DO  702  I *1  *  MM 
N«l*8 

702  PCOXl I )«PCOZ(N) 

DO  707  KW*l,MM 
EXi»MM-KW 

r  x  < KW >3pCOX ( KW) *YND(NX  )**EXl 

707  XND(NX)eXND<NX)*FX(KW) 

GO  TO  672 

708  DO  37  KL*l*MM 
EX1*MM“KL+1 

FX(KL>«PC0Z(KU*YND(NX)**6X1 

37  XND(NX)*XND<NX)*FX(KL) 

672  CONTINUE 

DO  95  KC»1.NNT 
I F ( XND ( KC ) ) 33  *  35  *  35 

33  IF(KC-1)38.38.51 

38  XND ( KC ) *  0 • 0 
GO  TO  35 

51  XND<KC)*XNU<KC-.1> 

35  IF(XND(KC)-1, 0)95*34,34 

34  XND  <  KC ) ■ *  9995 

95  CONTINUE 

FIND  AVERAGE  VELOCITY 

96  SUMV*0  « 

DO  114  L«1*NNT 
VL ( L ) *  V I *XND( L ) 

114  SUMV«SUM V+VL ( L ) 

VLAV«8UMV/UN 

PL1FT« (2 .3769E-3/2. )  *VLAV**2*3, 14167* { 3 ,/8 , ) *RP ( 1 1 ) **2/144 , 
PL  I  FT - -PL  I F  T 

V0L*4./3.*3.14159*RP( i J )**3 

DO  730  MU«1*4 

GO  TQ<731, 732*733*734), MO 

731  DNS«75. 

GO  TO  735 

732  DNS*90  * 

GO  TO  735 
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73$  dn$»1Q5, 

SO  70  735 
734  DNS® 125 » 

/«  PWf*DN S/1728, *VOL 
•  PRINT  735#  DNS 

7#8: FORMAT  <  /ABHPART.I QL§  &INI!  TY  <1>  .  . ~  ■ 

PLTOT®  PUirT-PWTtFl.SC 
PMAS®PWT/32.16 

plac.pl  tot/pmas  . . .  . 

!F<PLAC)84 2, 843.843 

842  PL AC*0  .  0 

843  PRINT  110  ,  , 

110  FORMAT < /22HAV6R A Q£  LOCAL  VELOCITY. 3X.16HUFT  ON  PARTICLE. 3X. 

1  26HN6T  LIFT  FORCE  ON  PARTI CLE, *X, 21HPART JCU:  ACCELERATION) 

PRINT  120 

120  FORMAT  <  7X.8HC FT/SEC) » 7x»22hDUE  TO  WALL  CONST RAlNT»3Qx.8lHlN  VgRTIC 
1AL  DIRECT  ION/ ) 

PRINT  840.VLAV.PLIFT.PLTQT.PLAC 
84q  FORMAT (F20.1O.E2U.  10. 5X.F20.*0.5X.F2|).  10/) 

730  CONTINUE 
GO  TO  115 

666  PRINT  667 

667  FORMAT (/5HERR0R) 

GO  TO  668  • 

115  CONTINUE 

668  RETURN  * 

END 

SUBROUTINE  HOOT ( NNT, MM. PCOZ. YND. XND.NCC > 

DIMENSION  YND(70).XN3(70).PCOZ(30).DVX(30>.FX(30) 

DO  ?10  L®1.NNT 
IF(NCC-24>199, 199,201 
199  IF(YND<L)-. 425)200, 201, 2Q1 
201  XND( L) ®0 . 1  ■  • 

106  SFX® 0 . 

SD VX ®0 « 

QOlOONal.MM 

MX1®MM-N*t 

MX2«MM"N 

EX1«MX1 

FX(N)«PCOZ(N)*XND<L)**MXi 
IF(N“MM)101, 102,101 

101  DVX(N)®EXl*PC02(N)*XND<L)**MX2 
GO  TO  103 

102  DVX(N)®PC02(N) 

103  SFX® SFX  +  FX ( N ) 

100  SDVX®SDVX+DVX<N> 

RMX®<SFX-YND(L) J/SDVX 

IF  < ABSF (KMX )- . 01  )  108 , 10m , 10  4 

104  XND<L)>XN0(L)-RMX 
MLM®1000/MM 

-  IF(XND<D*2***MLM)106. 204,204 

QO.  TO  210 
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108  IF(NCC-24>109#109,210 

109  IF(NNT-2)210,210,110 

110  IF (XND(L)~l>0)210»203i203 
253  IFJYNSXD-a.O  >204*204; 205 
2 0 5  XND ( U ) » . 9999 

GO  TO  210 

"8  o'*  jFi'L-1  V3i)4/364  3  0  9 

304  XND ( L ) *0  #  0 
GO  TO  210 

305  XND<L)»XND(L-1) 

GO  TO  210 

200  GO  TO  210 
210  CONTINUE 
RETURN 
END 
END 


DRAG  FORCE: 

Input . 

Program  Is  generated  in  the  same  manner  as  the  lift 
program. 

Intermediate  Results  and  Output. 

DELT  Boundary  layer  thickness 

XDEL  at  top  of  boundary  layer 

FDT  Drag  due  to  pressure  force 

FFDT  Drag  due  to  frloiton 

PWT  Parti ole  weight 

PDAC  Horizontal  aooeleratio$  of  partiole 


SUBROUTINE  DR AG < PCOZ , MM , NCC *  V J , RG , 2G , ON. YHVX , VMX ) 

DIMENSION  DP(10).UN(40)iRP(10)#rGU(40)#YGBC40)»YL(70). YG<70) 

1  AR(70).YND(70),XND(70).VL(70)#PCOZ<30).XH4o>»FX(30)* 

2  XDEL<5)#YMAX<5)#PCOX<10)iARF(7Q) 

OP l 1 ) 8 . 06 

OP 1 2  >  8  *  04 
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u  r  \  u  /  -  ■«.  u  £ 

DP(4)» ,01 
DP ( 5 ) 8  «  0  0  7 
DP(6)«.12S 

DP(7)»  ,  003  !' 

NN«3l 

RQN«RQ/DN 

ZGN8ZG/DN 

PBINT  All,  ON 

611  FORMAT  I22HN0ZZLE  DIAMETER  (FT)  ■  F5.2) 

MINT  612*  VJ 

“^it-f-tjrwt-^  -Ffvi'-'T-  ■■ 

PRINT  400* RON 

400  FORMAT (26HN0NDJMEN  RADIAL  DISTANCE  ■  F20.10) 

PRINT  401*  ZQN 

4d  FORMAT(28HNONDIMEN  HEIQHTH  OF.  NOZZLE  *  F2Q.10) 

00  115  1 1  «1 « 7 
PRINT  403*  CP(tl) 

4o3  FORMAT (  23HPARTICLE  DIAMETER-IN.  ■  F2o.1U  > 

RP< 1 1 ) "DPI  1 1 )/2 . 

UN  I *NN 

TK«RP( II )/<UNI-,5> 

YLC8Q. 

SUMAR-0. 

NATsNN-1 
DO  3  N«1»NAT  ■ 

YL(N)sYLC*TK/2, 

AM  =  N 

XL(N)*SORTF  (RPU  I?  *2-VL  <  N)  **2  ) 
AR(N)»2,*XL(N)*TX 
ARF<N)s4,*XL(N)*TK 
IF<N-1)6.5,8 

5  SUMAR»SUMAR^3. 1416/2, *ARF(N> 

GO  TO  6 

O  SUMAR*SUMAR+3.1416*ARF(N) 

6  YGU<N)»RP< I  I ) ♦YLC 
YQBC  N) »RP ( 1  I > -VLC 

3  YLC8 AM*TK 

QNU». 0000003745/. 0023769 
IF(NCC "25)302*303,302 
303  VI«VMX 

QO  TO  700 

302  IF(ZQ/DN". 5)640*640*661 
661  KZBSZG/DN 

QO  T0<640»64l, 643*643), KZQ 

640  IF(RGN". 35)60,60, 61 
61  IF(RGN-. 75)16*18,19 
19  IFIRQN-l. 2)20*20, 22 
60  VI8VJ*S3RTF(RGN**2*.Q175#RGN) 

QO  TO  700 

16  VI8VJ*SQKTF<-3,2*RGN**2*5.13*RGN-1.26) 

GO  TO  700 

2q  VI*Vj*SOKTF (-2,6B*RGN**2^5,46*RGN-1.82) 

GO  TO  700 

22  VIsVJ*SQRTF(.256*RGN**2-l146*RGN+2,3ll' 

GO  TO  700 

641  IF(RGN", 35)6o,6o, 644 

644  *  I F  < RGN-1. 0  )645, 645,647 

645  VIaVJ*SORTF(-,0985*RQN**2+l,l3*RGN".241) 

GO  TO  700 

647  VI«VJ*SQRTF(.077*RGN**2",641#RGN+1.354) 

GO  TO  700 
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64j  ImRuN-. 40)62, 62, 654 

654  IF(RGN-1. 2)655, 655, 657 

655  VI«Vj*SQHTFC“.5*RGN**2«,1.2*RGN-.3> 

GO  TO  700 

657  V I "V J*SQRTF (  .  q4o4*RGN**2- , 3l8*RGN* . 74* ) 
GO  TO  700 

62  V  I  • V J*SQHTF  (  ,625*Rgn**2«-.o25*RQN> 

700  NNT«2*NN-3 

CM'S  fUI«4  titi  T 

I F < N*NN ) 24 , 26* 26 

, .  14  JfMNiJ..YQB.(  NNrN.) .  .. 

. GO  TO  27  .  . 

26  YQ ( N ) ■YGU( N"NN*2 ) 

1 F(N-NN)27»27#28 
27  AR( N) * AR ( NN*N) 

GO  TO  977 

28  AR(N)«AR(N-NN*2) 

977  IF(NCC-26)99, 900,900 
900  AAA»VI/RG 

YND(N)«YG(N)*SQRTF(2.*AAA/GNU)/12. 

GO  TO  29  ' 

99  IF(NCC-25)978, 301,900 
301  YND(N)byq(N)/YHVX 
GO  TO  29 

978  YND(N) s<  YG<  N )/(12 .*ON) ) *SQRTF ( V J*DN/GNU) 
29  CONTINUE 
672  DO  805  NW=l,NNT 
XND ( NW ) s  o , 0 
IF(NCC-25)808,B04,808 
804  IF(YND(NW)-, 15)800, 800, 801 

800  MMbB 

DO  811  1 2*1 , MM 
811  PCQXt IZ)sPcOZ(IZ) 

GO  TO  808 
801  MM*  4 

DO  802  I  si , km 
802  PCQXCI )bPCOZ( 1*8) 

806  DO  807 
EX1*MM-KW 

FX<KW)bPCOX(KW)*YND(NW)**EX1 
807  XND(NM)sXND(NW)4>FX(KN) 

GO  TO  805 

808  DO  809  K Y«i *MM 
EXl»MM-KY+l 

FX(KY) «PC02(KY)*YND<NW)**EXl 
809  XND(NW)sxnD(NM)^FX(KY) 

805  CONTINUE 
75  DO  32  KC»i,NNT 

I F ( XND (KC))33»670,670 
33  1F<KC-1)72, 72,73 
72  XND ( KC ) * 0 i 0 
GO  TO  670 

73  XND(KC)«XND(KC-1) 

670  I F ( XND  « KC  > -1 • 0)32,32,34 
34  XNDCKC)*. 9999 
32  CONTINUE 

IF( I I"l)  74, 74 ,80 
74  IF ( NCC-25 ) 78, 76, 79 
79  DEUTb4, 6*1?./S0RTF(2, *AAA/GNU) 

GO  TO  80 
76  XDELsO.99 

CALL  ROOT  <1,MM,PC0X,XD6L, YMAX.NCC) 


GO  TO  80 
78  XDEL  < 1 ) » . 99 

CALL  ROOT (l,MM,PCO/»XD£L, YMAX.NCC) 
0ELTsYMAX<l>*DN/SQRrF<VJ*DN/aNU>*12, 

80  GO  TO  31 
31  D035Ka 1#  NNT 

35  VL(K)«VI*XND(K)  '  . 

82  NMID«(NN1 *1 )/? 

RN«VL(NM10)*DP(II)/<12.*ONU> 

I FIRN^lOO • ) 730  >730 # 731 
7  3 1  1 F  <  R  N  -10  0  0  0  <"1-7-32  ;  732*7  33 
730  CD*30./<RN**„699) 

GO  TO  734 

732  CDal.44/(HM*t>, 0396) 

GO  TO  734 

733  CD*1 • 18 

734  FDT*0. 

D036K«1,NNT 

FDa2.3769E»3*CD*VL<K)**2*AR<K)/2B8, 

36  FDTbFDT+FD 
ARS»4.*3.14159*RP< 1  I )**2 
ARDsAHS-SUHA'R 

AMD*  ARD/ AHS 
ACORsl ,♦ ARD 
CF *  .  664/SURTF ( RN ) 

FFDT  »  0 . 

0O40M«l»NNT 

FFD»2.3769E-3*CF*VL(M)**2*ARF(M)/268. 

40  FFDT *FF0  f +  FFD 
CFB«4.*CF 

TOL«ABSF(CKB*. 00001)  • 

CALL  CUBE(CF8.T0L»CFB3) 

CBDa  0 • 1/CFB3 
FBDT«0, 

DO  42  MD=1.NNT 

FBD«2.3769fc-3*CBD*VL<MU>**2*AR(MD>/2B8, 

42  FBUTaFBDT+FBD 
PRINT  38 

38  FORMAT  </4X  1 16HPRESSURE  DRAG"LB*4X»16HFRlCTIr‘N  DRAG-LB .  8*  ,  12HBASE  D 
1RAG-LB,2X,24hBOUNPARY  LAYER  TH  I CKNE5S, 2X , 31HSPHERE  APPROXIMATION  D 
2IFFERENCE//) 

PRINT  41.  FDT, FFDT, FBDT.DRLT, ARD 

41  FORMAT  <  3E20.10.2F20 ,10) 

V0L«4./3.*3.14159*RP< ll)**3 
DO  720  MO»l » 4 

GO  TO<721. 722.723, 736), MO 

721  DNS»75. 

GO  TO  724 

722  DNS«90 , 

GO  TO  724 

723  QNS*105 , 

GO  TO  724 

736  DNS*125. 

724  Pwl«DN S/1728, *VOL 
PRINT  725, DNS 

725  F0RMAT(/18HPARTICLE  DENSITY  «  F6.1> 

TOIL«(FOT*FFDT-,707*PWT)*ACOR 
PRINT  620,  TOTl 

620  F0RMAT(/34HNET  FORCE  IN  DIRECTION  OF  MOTION  ■  E20.10) 

PMAS-PWT/32.16 

pdac«totl/pmas 
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r 

IF(PDAC>842/843,843 

n  4  n  m  i%  4  0*  -  m  + 

Qffi  ruMU- V • U 

043  PRINT  821,  PQAC 

021  F08M11US6HPARTICLS  ACCELBRAT1 ON  IN  HORIZONTAL  PlRBCTIQN  <FT/fBfl) 
“  l'i  F2B.I0) 

720  CONTINUE 
115  CONTINUE 

IF ( NCC-24 ) 668# 668» 599 
899  CONTINUE 
660  RETURN 

END  . 
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